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Summary 

Models that predict the failure probability of monolithic glass and ceramic components under 
multiaxial loading have been developed by authors such as Batdorf, Evans, and Matsuo. These “unit- 
sphere” failure models assume that the strength-controlling flaws are randomly oriented, noninteracting 
planar microcracks of specified geometry but of variable size. This report develops a formulation to 
describe the probability density distribution of the orientation of critical-strength-controlling flaws that 
results from an applied load. This distribution is a function of the multiaxial stress state, the shear 
sensitivity of the flaws, the Weibull modulus, and the strength anisotropy. Examples are provided 
showing the predicted response on the unit sphere for various stress states for isotropic and transversely 
isotropic (anisotropic) materials — including the most probable orientation of critical flaws for offset 
uniaxial loads with strength anisotropy. The author anticipates that this information could be used to 
determine anisotropic stiffness degradation or anisotropic damage evolution for individual brittle (or 
quasi-brittle) composite-material constituents within finite-element or micromechanics-based software. 

1.0 Introduction 

The term “unit sphere” as used herein refers to the models that were developed by Batdorf and Crose 
(1974), Batdorf and Heinisch (1978), Evans (1978), and Matsuo (1981) to predict the probability of 
failure of brittle materials under multiaxial loading. These models use a unit radius sphere representing 
the random orientation of flaws to calculate the effect of multiaxial stresses on material reliability. This 
approach assumes that the strength-controlling flaws are randomly oriented, noninteracting planar micro- 
cracks of specified geometry but of variable size. Fracture mechanics relationships for mixed-mode 
(modes I, II, and III) crack growth — combined with the weakest-link theory and integration over the 
surface area of a unit radius sphere representing all possible orientations of microcracks — are used to 
calculate the material probability of failure. 

This unit-sphere methodology was originally introduced within a statistical theory of brittle material 
strength by Weibull (1939), though without consideration for the mechanics of crack growth. Nemeth 
(2013) extended this methodology to predict anisotropic (transversely isotropic) brittle material strength 
response and demonstrated how this could be used as a failure criterion inside finite-element (or alter- 
natively, micromechanics-based) software codes for individual brittle (or quasi-brittle) composite- 
material constituents to predict the overall strength response. 

The unique feature of the unit-sphere methodology is the assumption that the presence of flaws in a 
brittle material drives the failure response and that these flaws are microcracks with an assumed orien- 
tation and geometry. In this manner, the stochastic strength response of the material can be predicted for 
an arbitrary multiaxial stress state. An interesting consequence of this methodology is that the orientation 
of the critical flaws (the failure-initiating flaws) relative to the applied multiaxial stress state can be 
predicted as a probabilistic distribution. 

This report develops the generalized formulation to describe the probability density distribution 
(PDD) of the orientation of critical-strength-controlling flaws for the unit-sphere model. The derivation 
starts from the most basic level with the description of the weakest link mechanism and the Weibull 
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distribution. Following that is the description of the Batdorf unit-sphere methodology for the isotropic 
strength response. The extensions for the unit-sphere methodology with the anisotropic (transversely 
isotropic) strength response are provided to complete the description of all relevant parameters. These 
extensions are for (1) flaw orientation anisotropy, where a preexisting microcrack has a higher likelihood 
of being oriented in one direction than another direction, and (2) critical strength or fracture toughness 
anisotropy, where the level of critical strength a lr or fracture toughness K lc for mode I crack propagation 
changes with the orientation to the microstructure (for a tensile mode of failure only). This leads to the 
generalized formulation for the PDD of the critical flaw orientation over the unit sphere. Examples are 
provided to demonstrate this methodology for isotropic and anisotropic materials for various stress states 
and off-axis loading. 

If one knows the orientation distribution of critical flaws (and the subsequent direction of crack propa- 
gation), the anisotropic stiffness degradation (the anisotropic elastic constants associated with the damaged 
material) can be determined. The author anticipates that this will be helpful for a follow-on phase of this 
effort, not described here, of enabling the unit-sphere failure criterion methodology to work with NASA’s 
micromechanics analysis code/generalized method of cells (MAC/ GMC) (Bednarcyk and Arnold (2002)), 
which is based on the GMC family of micromechanics theories, including doubly and triply periodic versions 
of the GMC (Aboudi (1995)) and the High-Fidelity Generalized Method of Cells (HFGMC) (Aboudi et al. 
(2003)). This incorporation will allow the full exercise of the unit-sphere methodology, including incremental 
time/load steps and fatigue analysis (as described in Nemeth et al. (2005)), to predict the durability (strength 
and lifetime) of composite laminates and woven composite structures. 

The unit-sphere methodology provides an improved mechanistic basis to the problem of predicting 
the strength response of an anisotropic material under multiaxial loading in comparison to polynomial 
interaction equation formulations such as Tsai-Wu, Tsai-Hill, and Hashin, among others. A unique 
physical attribute of this model is the ability to predict the critical-flaw-angle distribution for an arbitrary 
stress state. This distribution is derived and demonstrated in this report. 


2.0 Model Description 

2.1 Weibull Distribution 


Consider a stressed component containing many flaws (microcracks), and assume that failure is 
due to any number of independent and mutually exclusive mechanisms (links). Each link involves an 
infinitesimal probability of failure. Discretize the component into n incremental links. The probabil- 
ity of survival (P s v), of the i th link is related to the probability of failure (P/v)i of the z th link by 
(P s v)i = [1 - (Pjv)i\, and the resultant probability of survival of the whole structure is the product of 
the individual probabilities of survival: 



where the subscript V denotes volume-dependent terms. Equation (1) arises from the approximation to the 
exponential series 

2 3 4 

XXX 

e x = 1 + x H 1 1 1 - ... (for all real values of x) (2) 

2! 3! 4! 


where {Pjv)i is assumed to be a small number such that the higher order terms in Equation (2) can be 
neglected. 

Assume the existence of a function q ( {a), referred to as the crack-density function, representing the 
number of flaws per unit volume that have a strength equal to or less than a. Under a local tensile stress 
a,, the probability of failure of the z th link, representing the incremental volume AT,, is 

(P/v)i=[r \vi<ydAV,] (3) 
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where the incremental volume AF,- is arbitrarily small such that the value of the expression within the 
brackets is much less than one. If one applies a uniform tensile stress g, such that a = a, for all 
incremental volumes A V,-, then from Equation (1) the resultant probability of survival for material 
volume V, where V is the sum of all AK„ is 


P sV = exp[-q,{a)E] 


( 4 ) 


Equation (4) can also be derived from the Poisson probability density function. This distribution is 
described by (see Hoel et al. (1971), for example) 


X x exp(-k) 


P(X = x) = f{x) = \ 


x\ 


[0 


x = 0,1,2, ... 

elsewhere 


(5) 


where X is a positive number. The real- valued function f(x) is the discrete density function of random 
variable X where P(X = x ) is the probability that a discrete real- valued random variable X equals a 
possible value x. The Poisson distribution approximates the binomial distribution for large values of n, 
where n is the number of Bernoulli trials with success probability p = Xln at each trial. Equation (4) is 
obtained when P(X= 0) is computed for n = V and p = rp{a); hence, 

P S V = p ( x = 0) = ^ eX j! ( = exp(-X) = exp[- %(g) v] (6) 


Equation (6) calculates the probability that no flaws of strength a or less are present in the material 
volume V and, therefore, represents the survival probability of the material under applied load a. 

The probability of failure for the uniformly stressed volume V is 

P/v =l-exp[-r| F (a)E] (7) 


where V is the total volume. 

Weibull introduced a three-parameter power function for the crack-density function q ( (a): 




( \ m v 

V-VuV 

a oV 


( 8 ) 


where G^is the threshold stress parameter. The scale parameter a 0 y then corresponds to the stress level 
where 63.21 percent of tensile specimens with unit volumes would fracture. The scale parameter has 
dimensions of stress x (volume) Umy , where m v is the shape parameter (Weibull modulus) — a dimension- 
less parameter that measures the degree of strength variability. As m v increases, the dispersion decreases. 
The threshold stress parameter a u y is usually taken as zero for brittle materials (ceramics and glasses). 
This parameter is the value of the applied stress below which the failure probability is zero. When this 
parameter is zero, the two-parameter Weibull model is obtained in Equation (7). The two-parameter 
crack-density function is expressed as 


%(c0 = 


( \ m v 

G 


V J 


= k 


wV G 


( 9 ) 


where k w v= (a 0 v) mV is the uniaxial Weibull crack-density coefficient (for an applied uniaxial load). 
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2.2 Batdorf Unit-Sphere Model for Isotropic Strength Response 

Batdorf and Crose (1974) propose a statistical theory in which attention is focused on cracks and their 
failure under stress. Flaws (assumed to be microcracks in shape and size) are taken to be uniformly distri- 
buted and randomly oriented in the material bulk. Fracture is assumed to depend only on the tensile stress 
acting normal to the crack plane; hence, shear insensitivity is inherent to this first model. Subsequently, 
Batdorf and Heinisch (1978) included the detrimental effects of shear traction on a flaw plane. Their 
method applies fracture mechanics concepts by combining the crack geometry and a mixed-mode fracture 
criterion to describe the condition for crack growth. 

In the Batdorf unit-sphere theory, the incremental failure probability AP fV under an applied multiaxial 
state of stress X at a given location in the component can be described as the product of two probabilities: 

APfvfcvieqoWh&PwPlV (10) 

where A P lv is the probability of the existence in incremental volume AFof a crack having an equivalent, 
or effective, critical stress between o leqc and ai eqc + Aa leqc . Critical stress a k is defined as the remote, 
uniaxial fracture strength of a material containing a given crack in mode I loading. The term a leqc denotes 
an effective (or equivalent) critical mode I stress from applied multiaxial stresses. The second probability, 
P 2 v, denotes the probability that a crack of critical stress G\ eqc will be oriented in a direction such that an 
effective stress <3i eq (which is a function of fracture criterion, stress state, and crack configuration) satisfies 
the condition G\ eq > G\ eqc . The effective stress a ieq is defined as the equivalent mode I stress that a flaw 
would experience when subjected to a multiaxial stress state that results in mode I, II, and III crack sur- 
face displacements, and a leqc is the threshold value of a Te? where unstable catastrophic crack growth 
ensues. An incremental volume AK is used in Equation (10) because an infinitesimal volume d V cannot 
enclose a crack of critical stress 0 \ eqc and associated critical crack length a c . 

The effective stress 0 ] eq represents an equivalent normal stress on the crack face from the combined 
action of the normal stress ct„ and the shear stress x. The microcrack orientation is defined by the angular 
coordinates a and [1, where the direction normal to the plane of the microcrack is specified by the radial 
line defined by a and p in space (see Fig. 1(a)). For the sake of brevity, the development of the effective 
stress equations is not shown (for details, see Nemeth et al. (2003, 2005)). For a tensile mode of failure 
assuming a penny-shaped crack with the Shetty mixed-mode fracture criterion (Shetty (1987)), the 
effective stress becomes 




-| 

2 

1 


4x 


Gleq ~T 

2 

o„+,c>Z + 

T(2-v)J 



( 11 ) 


where v is Poisson’s ratio and C is the Shetty shear-sensitivity coefficient, with values typically in the 
range 0.80 < C < 2.0. As C increases, the response becomes progressively more shear insensitive. Shear 
increases the equivalent stress as shown in Equation (11), and this has a deleterious effect on the predicted 
material strength. For a penny-shaped crack with a material having a Poisson’s ratio v ~ 0.22 and C = 0.80, 
0.85, 1.05, and 1.10, Equation (1 1) approximates, respectively, the following criteria: Ichikawa’s maximum 
energy-release-rate approximation (Ichikawa (1991)), the maximum tangential stress (Erdogan and Sih 
(1963)), the maximum strain-energy-release-rate formulation of Hellen and Blackburn (1975), and colinear 
crack extension. The value of C also can be fit empirically to experimental data — either on introduced 
cracks (as is done in Shetty (1987)) or on specimens being tested multiaxially. 
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Figure 1 . — Projection of equivalent stress onto a unit radius sphere. 

(a) Cauchy stress components on an infinitesimal tetrahedron 
resolving a normal stress, a n , and a resultant shear stress, x, on a 
plane normal to the direction defined by angular coordinates a and p. 

(b) Projection of equivalent stress onto a unit radius sphere in the 
global coordinate system. The unit radius sphere represents all 
possible flaw orientations, where ai eq is an equivalent (or effective) 
stress; a x , a y , and a z are normal orthogonal stress components; and 
ixy, Tyz, and izx are shear stress components. An infinitesimal area, dA, 
on the unit sphere represents a particular flaw orientation (a direction 
normal to the flaw plane), and c\ eq is a function of a„ and i for an 
assumed crack shape and multiaxial fracture criterion. 
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A material is assumed to have n individual failure modes. Material reliability (probability of survival) 
is assumed to be the product of the survival probability of all of the failure modes. These failure modes 
can be from different flaw populations, from the flaws associated with different material constituents 
(i.e., composites as described in Nemeth (2013)), or from different physical mechanisms that govern the 
material failure such as tensile loading versus compressive loading. Here the compressive mode of failure 
is assumed to be of a different nature than the tensile failure mode — possibly involving the interaction of 
arrested cracks (flaws with initial crack growth that subsequently stops because of the angle of the crack 
growth and the interaction with the local stress field). Regardless, the Weibull distribution here is 
assumed to describe the stochastic strength response phenomenologically. This is argued further in 
Nemeth and Bratton (2011). The compressive failure mode is assumed to be controlled by shear stress, 
and a simple Tresca-like effective stress o Ie? relation can be prescribed as 


Gleq = 2x (12) 

The multiplier of 2 in Equation (12) was chosen so that the maximum effective stress on the unit sphere in 
pure uniaxial compression is equal to the applied compressive stress. When the normal stress component g„ 
on the crack face is tensile, then the value of <3\ eq in Equation (12) is set to zero in the unit sphere evalu- 
ation. A more rigorous treatment to account for compression would include the frictional effects of the 
opposed crack surfaces in contact in the effective stress relation. However, Equation (12) is sufficient for 
the demonstration purposes of this report. 

2.2.1 Relationship of Crack Probability of Occurrence AP\ v to the Crack-Density Function r\y 

The strength of a component containing a flaw population is related to the critical flaw size, which is 
used implicitly in statistical fracture theories. Batdorf and Crose (1974) describe the crack probability of 
occurrence AP ]V as 


AP W =AV 


dr lp(gi ggc ) 


do 


I eqc 


leqc 


(13) 


This is better understood using Figure 2. In Figure 2(a), for example, a point on the curve r\t(o\ eqc )AV for 
the value Oi e?c j represents the number of flaws in AK having a strength equal to or less than <Ji eqc ,i. Then, 
from Figure 2(a), 

T\vi?leqc,2)kV - q ,{G I(?9 C ,i)AK 

= Probability that a flaw of strength o u , q(q2 or less is in A K 

- Probability that a flaw of strength a re?c j or less is in A K (14) 

= Probability that a flaw of strength between Oi egc ,i and 0\ eqCt 2 is in AK 

Equivalently in Figure 2(b), 

[r\A<5\eqc) + dr\y]AV- r\f{aieqc)AV 

= Probability that a flaw of strength (<5\ eqc + doi ( , (/( ) or less is in AK 

- Probability that a flaw of strength <3\ eqc or less is in AK (15) 

= Probability that a flaw of strength between a leqc and (o k?c + dO| ( , (/( ) is in AK 

= AK dr| , 


The slope s of the relation r|i{G|„ /( ) AK is s = (drp /dai e?c ) AK as shown in Figure 2(b); therefore, 
multiplying the slope s by doi eqc yields the result of Equation (13). 
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Figure 2. — Crack-density function, r|i/, versus the critical equivalent 
stress, oieqc, where A\/is incremental volume, shown to explain 
Equation (13). (a) Incremental differences in values, (b) The same 
differences described in terms of infinitesimal values. 


2.2.2 Random Flaw Orientation as Described by a Unit Sphere 

The term P 2V shown in Equation (10) is expressed as 


P 2V - 



(16) 


where Q(Z, <jy qc ) is the area of the solid angle projected onto the unit radius sphere in stress space (see 
Fig. 1(b)) containing all the crack orientations for which a leq > <5{ eqc for the applied far-field multiaxial 
stress state I. The infinitesimal area, cL4, on the unit sphere represents a particular flaw orientation 
(a direction normal to the flaw plane), and e> u , q is an equivalent stress, which is a function of an assumed 
crack shape and multiaxial fracture criterion. The constant An is the surface area of a unit radius sphere 
and corresponds to a solid angle containing all possible flaw orientations. Equation (16) assumes that the 
flaws are randomly oriented with a uniform distribution (equal likelihood of being oriented at any angles 
a and [1). If it is assumed that the other properties like crack geometry and K lc are also independent of 
a and [1, the strength response becomes isotropic (independent of the applied direction of load). 

2.2.3 Unit-Sphere Failure Probability Formulation for Isotropy 

Integrating Equation (10) with respect to Gi eqc and using Equation (4) yields the probability of 
survival of a (uniformly stressed) volume element AT,: 
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(17) 


); = expj-AF 

where is the maximum effective stress that a randomly oriented flaw could experience from the 

given stress state. Hence, integrating over the entire volume of the component results in the failure 
probability: 



f ®leq,i 

'o 


®Ieqc ) d fl V ^Pleqc ) 


47T 


d ®Ieqc 


d ®leqc 




fl-exp -J, 


) d fl V ( CT I eqc ) 


4 n 


' eqc 


eqc 


dF 


(18) 


where G[„ /-max is the maximum effective stress that a randomly oriented flaw could experience in the 
component and X is allowed to vary within the component. 

Because of the effect of multiaxial stress states, the crack-density function rp/(ai e?c ) is altered in the 
Batdorf formulation. Batdorf and Heinisch, (1978) have typically used a power-law formulation for 
flr(Gi eqc), which leads to the Batdorf crack-density function of the form 


% \V\eqc ) = k BV^leqc ( x > z > a > P) : 


^BV^leqc ^’ y,Z ’ a, P) 


>oV 


(19) 


where x, y, and z correspond to the location within the component and a and P are the unit-sphere 
orientation angles. The term k BV is the Batdorf crack-density coefficient, and k BV is the normalized 
Batdorf crack-density coefficient whereby k BV = k BV /k wV as described in Nemeth et al. (2003, 2005). 
The material Batdorf crack-density coefficient k HV and the Weibull modulus m v are evaluated from 
experimental inert strength fracture data. By convention, k BV is calibrated to a uniaxial stress state. 
Batdorf and Crose (1974) initially proposed a Taylor series expansion for r|i(G| <a/a ), but this method has 
computational difficulties. Note that p i(oi eqc ) has units of inverse volume. 

Although the Weibull (Eq. (9)) and Batdorf (Eq. (19)) crack-density functions are similar in form, 
they are not the same. The Weibull function simply depends on the applied uniaxial stress distribution g 
and is the only term other than m v and the volume necessary to calculate P,v. It does not extrapolate to 
other multiaxial stress states. The Batdorf function depends on the equivalent mode I strength of the 
crack a leqc , which is probabilistic and must be integrated over a range of values for a given stress state. 
Furthermore, to obtain P fV , a crack orientation function, P 2 v, must be considered in addition to the density 
function and the volume. Finally, the Batdorf coefficient k B y cannot be calculated from inert strength data 
until a fracture criterion and crack shape are chosen — in contrast to the Weibull coefficient k w y, which 
depends only on normalizing or calibrating to experimental rupture data. 

To determine a component probability of failure, one must evaluate P 2 v (Eq. (16)) for each elemental 
volume A V„ within which a uniform multiaxial stress state I is assumed. The solid angle Q(I, G la/C ) 
depends on the selected fracture criterion, the crack configuration, and the applied stress state. For 
multiaxial stress states, with few exceptions, Q(X, 0 \ eqc ) must be determined numerically. For a sphere of 
unit radius (see Fig. 1(b)), an elemental surface area of the sphere is cL4 = sin a dp da. If we project onto 
the spherical surface the equivalent (effective) stress Gi ( , ? (X, a, P), the solid angle Q(X, o ]a/c ) will be the 
area of the sphere containing all the projected equivalent stresses satisfying a ]eq > a leqc . Note that the 
symmetry of a u , q means that only one-eighth of the unit sphere needs to be integrated for an isotropic 
strength response. However, for transverse isotropy, one-half of the sphere should be integrated. When 
one -half of the unit sphere is considered, 
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where 


^(^’ CT l£Y/c) f 1 ") f 2 71 f Jt/2 / \ . 

4n = 1 2 ^ J J 0 Jo H{o leq ,o leqc )smadadp 

H(C>i eq, Gleqc) 1 Gj eq — C [ e qc 

C| t i y. Gi eqc) 0 G\eq ^ Gleqc 


The Heaviside function H(a \ eq , G\ eqc ) is a function of fracture criterion, stress state, crack configura- 
tion, and crack orientation. Substituting into Equation (18) and integrating first with respect to G] eqc for 
each individual angle a or P on the unit sphere changes the component failure probability to (Batdorf 
(1978a, b)) 


Pjy = 1 - exp 


\ v \T\o 2 ^v ) sin a da d P d F 


(21) 


where 


n V (<H eq )= k BV°Z V q ( X ’ T’ Z ’ P) 


For a given incremental volume, Ci eg (x, y, z, a, P) is the projected equivalent stress over the unit radius 
sphere in coordinate stress space as shown in Figure 1. 

Equation (21) circumvents the involved numerical integration of Q(I, a leqc ) as was originally 
developed in Nemeth et al. (1990). Equations (18) and (21) are equivalent formulations; however, 
Equation (2 1 ) is more convenient for computational purposes with the few exceptions where special- 
case closed-form solutions exist for particular stress states, crack geometries, and fracture criteria. 

2,3 Unit-Sphere Model for Transversely Isotropic Strength Response 

Two different physical mechanisms were considered in order to extend the unit-sphere model to 
account for anisotropic strength response: (1) flaw orientation anisotropy, whereby a preexisting micro- 
crack has a higher likelihood of being oriented in one direction than another direction and (2) critical 
strength or fracture toughness anisotropy, whereby the level of critical strength g Jc or fracture toughness 
K lc for mode 1 crack propagation changes with regard to the orientation of the microstructure (for a tensile 
mode of failure only). Flaw orientation anisotropy was previously considered by Buch et al. (1977), and 
critical strength anisotropy was previously considered by Batdorf (1973). Both models were developed to 
simulate the graphite anisotropic strength response. Extensions to these models, described and demon- 
strated in Nemeth (2013) for transversely isotropic strength response, are reproduced herein because 
terms from these formulations are used in the flaw orientation probability density function (PDF) 
formulation. The extensions include shear sensitivity for flaws and an improved functional form for the 
anisotropy equations. 

2.3.1 Flaw Orientation Anisotropy 

Flaw orientation anisotropy refers to the situation where a flaw has a higher likelihood of being 
oriented in one direction than another for a given critical strength. This means that a material will be 
stronger on average in one direction than another. An isotropic brittle material is equally strong in any 
direction, and thus its flaws are uniformly randomly oriented. However, in components made by 
processes such as extrusion or hot pressing, which induce texture, a bias will exist in the orientation 
distribution of processing flaws. Also, components finished by surface grinding will contain machining 
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damage in the form of surface cracks that are oriented parallel and transverse to the grinding direction. 

In composite materials, the interface or interfacial layer between the fiber and the matrix can act as a flaw 
with an orientation bias that induces an anisotropic strength response. 

In the following unit-sphere formulation, consideration is given to a flaw that has a higher likelihood 
of being oriented in one direction than another. If an argument analogous to Equation (10) is used, the 
incremental probability of failure under the applied multiaxial state of stress I for a given a and P orien- 
tation of the (direction or vector) normal to the particular flaw plane can be written as the product of two 
probabilities: 


A Pfi^Z, a, p, 0 [ e g C , AV) = AP lv p(a, P) da dp H(a 1eq , o eqc ) (22) 

where A/ J , v is the probability that a crack having an equivalent critical stress between 0 \ eqc and 
a leqc + &<5\ eqc exists in an incremental volume element AV The second probability, p(a, P) da dp, 
denotes the probability that a crack of critical stress a \ eqc is oriented in the range between a and (a + da) 
and between P and (P + dp), where the orientation of the microcrack is described by the vector that is 
normal to the plane of the microcrack. The Heaviside function H( c leq , a leqc ) is identical to that as defined in 
Equation (20), where H(a leq , a leqc ) = 1 when CTi,,,/!, a, P) > <3\ eqc and where H(c leq , o leqc ) = 0 when 
ai e9 (X, a, P) < Gi eqc . The effective stress Gi c? (X, a, P) is defined as the equivalent mode 1 stress oriented at 
coordinates a and P that a flaw would experience when subjected to a multiaxial stress state X that results 
in mode 1, 11, and 111 crack surface displacements, and o leqc is the threshold value of a leq where unstable 
catastrophic crack growth ensues. 

Analogous to Equations (16) and (20) for the unit sphere, 

P 2V = J 0 (3(a,p)f/(a Ie? ,G Ic )sina da dp (23) 


where 


£>(a,P) 


C(a,P) 

C 2it f tc/2 , , 

j 0 J 0 C(a,P) sin a da dp 


(24) 


and 


G| i ( y. Gif,) 1 


*5\eq — Gj ( 


H(<5i eq , G I , ) 0 G| j( y < G | 


Equation (24) is explained further in Section 2.4. The function £( a, P) describes the degree of anisotropy 
of the flaw orientation, where the normal direction to the flaw plane is given by angles a and p. 
Equations (23) and (24) modify Equation (21) for probability of failure: 


Pfy= 1 - exp 


Un>P) eg) Sin a da dp dV 


(25) 


where 
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riV (°Ieq h k BV °Z'q ( X ’ Y» Z ’ a ’ P) = k BV 


°Z V q (x,y,Z,CX,p) 
a oV 


\ m V 


The normalized crack-density coefficient k BV now includes the effect of flaw orientation anisotropy and 
is evaluated numerically from Equation (25) for a uniaxial stress state applied in one of the material 
coordinate system axis directions. 

For a transversely isotropic strength response, C, in Equation (24) is only a function of a. Buch et al. 
(1977) introduced a cosine power function for C,(a) = [cos (a)] 11, where <|) is a constant. This relation was 
modified in Nemeth (2013) to enhance the functional flexibility and is expressed as 


C( a ) 


- 

( M 


an 

cos 



v 2 A l J ~ 


0<a<A z 


£ (a) = 0 


A l < a < - 
2 


(26) 


Alternatively, 


<;(«)= o 


C( a ) 


( | 


- 

/ 

. 

71 


n 



a- 

— 

l 1 

[2A z, 


U 



n 

0 ft c x < A t 

2 


for — -A r < a < — 
2 2 


(27) 


where A and <|) are constants (with subscripts L or 7) that control the degree of anisotropy, and 0 < A < ti/ 2 
and <t> > 0. When A = ti/ 2 and ()> = 0, the isotropic strength response is obtained. Equations (26) and (27) 
are defined for one-half of the unit sphere (the top half, as shown in Fig. 1(b) where 0 < a < ti/2). 
Referring to Figure 3, Equation (26) represents the “polar-cap” or longitudinal L distribution of flaws and 
Equation (27) represents an “equatorial-belt” or transverse T distribution of flaws. The polar-cap distri- 
bution describes crack planes that are distributed symmetrically (centered) about a plane (in this case, the 
CT y -<7 Z plane), and the equatorial-belt distribution describes crack planes that are distributed symmetrically 
(centered) along a line (in this case the g x axis). For the polar-cap distribution, k BV is evaluated for a 
uniaxial stress along the g x direction; and for the equatorial-belt distribution, it is evaluated in the a y 
(or optionally the a z ) direction. 

The separate polar-cap and equatorial-belt distributions are introduced to describe individual and 
distinctly different failure modes. For a unidirectional fiber-reinforced composite, the polar cap can be 
used to represent the fiber strength distribution and the equatorial belt can be used to represent the matrix- 
fiber interface. The equatorial-belt and the polar-cap distributions can be considered to equivalently 
represent global failure planes, which are referred to as “action planes” in the Puck multiaxial strength 
model for composites (see Lutz (2006) for a description). The angular width of the belt or cap with regard 
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Figure 3. — Unit sphere with probability distribution functions (PDFs) describing the anisotropy of flaw 
orientation, where a and p are angular coordinates and A L and A r are constants in the flaw-orientation 
anisotropy function representing longitudinal and transverse distributions. Orientation is described with the 
normal to the crack plane. In this figure, two orientation functions are described: (1) a polar-cap distribution 
describing crack planes distributed symmetrically (centered) about a plane (in this case the <jy-a z plane) 
and (2) an equatorial-belt distribution where crack planes are distributed symmetrically (centered) along a 
line (in this case, the a x axis). 
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to angle a indicates the maximum extent of scatter of a flaw orientation. Flaws or action planes oriented 
outside of the scatter bands of the cap and belt are assumed to not exist or not contribute to the likelihood 
of failure. An example of this methodology is provided in Nemeth (2013) for a unidirectional polymer matrix 
composite (PMC) under biaxial loading. 


2.3.2 Critical Strength or Stress Intensity Factor Anisotropy 


Batdorf (1973) approached strength anisotropy using the G k strength ellipsoid approach (describing 
an ellipsoid rather than a unit sphere). This was applied to graphite — a mildly anisotropic material. In this 
report, strength anisotropy is used with the critical stress-intensity factor (fracture toughness) K lc varying 
with orientation angle on the unit sphere. This is functionally equivalent to a k varying with orientation, 
and it is also functionally equivalent to the size of the flaw changing with the orientation angle. With 
regard to a ceramic matrix composite (CMC), where failure from loading in the fiber direction is by 
matrix cracking with large-scale fiber bridging, fracture toughness cannot be defined on the global scale 
of the structure. In that case, one has to consider critical strength as a metric. However, it is acceptable to 
use fracture toughness on the local scale at the crack tip, where micromechanics can account for the 
bridging explicitly. In this report, fiber bridging is not considered directly. Of first-order importance 
herein is that local fracture toughness could change with the orientation of the flaw plane (or action plane) 
and the applied loading. The specific micromechanics of how this might occur is not considered here. 

The modeling approach taken here is similar to that described previously for flaw-orientation ani- 
sotropy (Section 2.3.1). The critical strength a k is defined as the fracture strength of the crack in mode 1 


loading and is proportional to K lc . Therefore, for anisotropic 


Kiln. Rt = (c. n,Ja. Rrt = 






where c is a constant (c = Yy[a^ with crack-shape geometry factor Y and critical crack length a c ), a k , n 


is the maximum value of a k over the unit sphere (for all a and P), and f lc (a,P) is a normalized function 

expressing the degree of this anisotropy. The unit-sphere formulation described previously is used, except 
the Heaviside step function in Equations (20) and (23) is modified as follows: 




_ CT I eg > 

/lc(a,P)“ 


®leqc ,max 


eq 

/ic(wP) 


< ®leqc, max 


(28) 


where a, eqc is substituted for a k to indicate generalized mixed-mode loading. 

For a transversely isotropic response where anisotropy is only a function of angle a, the function 
/ Ic (a) over the top half of the unit sphere (0 < a < n/2) is arbitrarily defined for the L distribution as 


/lc(a) = l 


m i 

— H COS 

an ^ 

y l 

[■- 1 ! 

U J 2 



l r Lj 


/lc(a) =1 


0<a <^ L 


> 


t,L< a< 


71 

2 


(29) 


and for the T distribution 
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/ic(a)= 1 


( 30 ) 


0 < a < E T 

2 


/lc( a ) = 1 



n 

2 


E, T <a<- 
1 2 


where E,l, £,/, yz, yr, r r.,, and r T are constants. The L subscript relates to the polar-cap strength anisotropy 
distribution for crack planes symmetrically distributed (centered) about a plane (in this case, the a y -a z 
plane), and the T subscript relates to the equatorial-belt strength anisotropy distribution for crack planes 
symmetrically distributed (centered) along a line (in this case the g x axis) (see also Fig. 3 for a reference 
frame). Figure 4 shows a schematic of Equations (29) and (30) in stress-strength space. 


°lc,x 

t 



Figure 4. — Schematic of Equations (29) and (30) in normalized 
strength space. Equation (29) correlates to the polar-cap 
distribution, and Equation (30) correlates to the equatorial-belt 
distribution, as shown in Figure 3. Here, El, Et, y l, and yrare 
constants in the critical mode I stress intensity anisotropy 
function representing longitudinal and transverse distributions; 
/j c (a) is the normalized anisotropy function of the critical 
mode I stress-intensity factor, K| C , or critical mode I strength, 
CT| C , as a function of angle a; and ct !c x , 0 | C y , and ct !c z are the 
normalized (by ai c , max ) orthogonal critical strength components 
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2.3.3 Generalized Unit-Sphere Failure Probability Formulation for Transverse Isotropy 

Equations (26), (27), (29), and (30) are arbitrary functions chosen for their flexibility to fit to data. 
Other functions could be easily substituted. For an isotropic distribution of the orientation of the flaws, 
p( a, P) is 1/271, which is consistent with Equation (21). Although not detailed herein, the Weibull 
modulus can also be made anisotropic, varying with angles a and p, by using, for example, the 
normalizing Equations (29) and (30) to modify m (dividing m by f\ c (a)). Therefore, in its most 
generalized form for transverse isotropy (considering flaw and critical strength orientation anisotropy), 
failure probability can be expressed as 


Pjy = I - exp 


• r 2 tu r it/ 2 

fJo Jo Jo 


^\v\ a \eqc) 


da 


I eqc 


p(a,P) sin a d a Ie9C da dpdE 


(31) 


and, if a power law is assumed for the crack-density function, 



)-^SF( a >P) 


7 Ie(? (x,y,z,a,p) 

/ic(wP) 


m v { a,P) 


= k 


BV 


a Ie<? (x,y,z,a,p) 

/ Ic .(a,p)a oV (a,p) 


N m F (a,P) 


(32) 


where critical strength or K lc orientation anisotropy is accounted for with the addition of the term f\ c . The 
Weibull modulus m v is allowed to be a function of a and p. Note that the Weibull scale parameter a oV is 
used in the crack-density formulation and that it is also allowed to be a function of a and p, but that it is 
assumed to be constant with spatial location x, y, z. Equation (3 1 ) can then be more simply expressed as 


Pjy = 1 - exp 


u„ 2 t; /2 *a 


%Pj7,a,P) 

/t c (a,p)a oF (a,p) 


m v ( a.P) 


p(a, P) sin a da dp d V 


(33) 


2.4 Orientation Distribution of Critical Flaws Under Multiaxial Load 

The unique feature of the unit-sphere methodology is the assumption that the flaws that are inherently 
present in a brittle material will drive the failure response and that these flaws are microcracks with an 
assumed orientation and geometry. In this manner, the stochastic strength response of the material can be 
predicted for an arbitrary multiaxial stress state. An interesting consequence of this methodology is that the 
orientation of the critical flaws (the failure-initiating flaws) relative to the applied stresses can be predicted 
as a probabilistic distribution. The PDF of the orientation of the critical fracture-causing flaws can be 
obtained from Equations (10), (13), (16), and (19) for an isotropic strength material; from Equations (22), 
(23), and (24) for flaw orientation anisotropy; and from Equation (28) for critical strength (fracture- 
toughness) orientation anisotropy. How this is obtained is described subsequently for transverse isotropy. 

For a given flaw (a flaw that is assumed to exist at some orientation a and P) and considering all 
possible flaw orientations, the likelihood must be 1 (in other words, 100-percent certainty) that the flaw 
exists at some orientation a and p. In other words. 


Probability that a flaw exists at some a and P -l~cy 


j- n j o n/ C (a,P)sinadadp 


(34) 


where c z is a normalizing constant and the function C(a, P) describes the degree of anisotropy of flaw 
orientation at angular coordinates a and P (where the normal direction to the flaw plane is given by angles a 
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and (3). The function C was defined in Equations (26) and (27) for a transversely isotropic strength response 
(which is expressed only as a function of angle a). Therefore, in Equation (34) c z must equal 


cz 


1 

\ITJ 2 <^( a >P) sina da dp 


(35) 


The probability that the flaw is oriented in the range between a and (a + da) and between P and 
(P + dp) is given by p(a, P) da dp in Equation (24). In this case, an infinitesimal area element cL4 = da 
dp on the unit sphere is represented as shown in Figure 5. In contrast to Figure 1(b), the sine term is not 
used here to describe cL4. Figure 1(b) maintains a convention previously used in other references such as 
Nemeth et al. (2003, 2005). The sine term is only necessary when integration over the whole unit sphere 
is required. In order to determine the PDF of the orientation of the critical flaws, cL4 must be kept constant 
versus angular coordinates a and p. If dA = sin a dp da were used, the nonsensical condition would be 
obtained that the probability of a flaw being oriented at a = 0 would always be zero regardless of the 
stress state and the magnitude of the stress. 

Equation (24) is the PDF of the degree of anisotropy of the flaw orientation. Integration of the 
numerator of Equation (24) over the entire unit sphere for <14 = sin a dp da yields a value identical to the 
denominator, satisfying the condition for a PDF that integration over all possible values must equal 1 (for 
the cumulative distribution function). 


«x 



Figure 5. — Infinitesimal area element, d4, on the unit sphere. 
In this case, d4 = da dp for a given a and p orientation of 
the normal to the flaw plane, where a\ eq is the effective 
stress. 


The probability that a flaw of strength <3\ eqc or less that is within the incremental volume element AV, 
where ai e9 (X, a, P) > G\ eqc , and that is oriented between a and (a + da) and between p and (P + dp), will 
fail is (see also Eq. (3)) 
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PfY&, a,p,AF) = p(a,p) dadpr| F (a,p) AV 


(36) 


[^(a,p) da dp ri (/ (a,P) AV 

r 2 k r Jt/2 

j 0 J 0 C(a,P)sma da dp 


where / J /j{X, a, P, AF) is expressed for a constant infinitesimal area element cL4 = da dp on the unit 
sphere as shown in Figure 5 regardless of the specific values of angles a and p. 

The failure probability for the whole unit sphere for incremental volume element A V is 


PJVV, AV): 


<;(a,P)r|,/(a,P)si n adadp AV 

Jo 2 X /2 ^( a ’P) sina da dp 


(37) 


Reexpressing Equation (37) using the crack-density function of Equation (32), using Equation (35) for 
c z , and considering critical strength or K lc orientation anisotropy yields 


Pjv(Z,AV) = 


So n So 2c zZM)k BV 


y leq 


(a,P) 


flc (a, P) a oF (a,p) 


m v ( a,P) 


sin a da dp 


AV 


(38) 


The probability that the critical flaw (the failure- initiating flaw) is oriented in the range between a 
and (a + da) and between P and (P + dp) that is under a given stress state X is found by dividing 
P t j/(X, a, P, AV) in Equation (36) by / J /( {X, AV) in Equation (37): 


, [C(a,P) da dp]r| F (a,P) 

p c (z.,a,p) , 2 

J J <^(a, P) q F (a, P) sin a da dp 


(39) 


The numerator in Equation (39) refers to an infinitesimal area element cL4 = da dp on the unit sphere as 
shown in Figure 5. The PDF of Equation (39) is therefore 


p z (X,a,P) = 


C(a,P)r|[/(a,P) 


r2n rn/2 

J q J q <;(a,P)r| F (a,P)smadadp 


(40) 


Integration of the numerator of Equation (40) over the entire unit sphere for cL4 = sin a dp da yields 
a value identical to the denominator, satisfying the condition for a PDF that integration over all possible 
values must equal 1 (for the cumulative distribution function). Substituting Equation (32) into 
Equation (40) yields 


£(a,P) k BV 


£? z (X,a,P) = - 


r . \«v(a,P) 

a l eg ( a ^) 

]\c ( a > P) 
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e2n rn/2 

Jo Jo ^ a ^ k BV 


f , n A%(a.P) 

a l eg ( a ^) 

/lc( a 5 P) 


(41) 


sin a da dp 
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However, for computational convenience, the effective stress a ]eq should be normalized by the scale 
parameter a 0 y and the normalized Batdorf crack-density coefficient k BV . If it is assumed that a oV and m v 
are independent of angles a and |3, 


C(a,P) 


^(2,a,P) = 


gi eg (a.P) 
/i c (a,P) c, oV 


n; /2 «-p) 


/ \m v 

CT feg( a >P) 

/ic(a,P) o oV 


(42) 


sin a da dp 


Equation (42) describes the PDF of the orientation of critical flaw normals from the applied stress 
state 2. However, these equations only describe the probability of the orientation of the critical flaw, they 
do not describe the angle of subsequent crack propagation. For noncoplanar crack growth, the action of a 
simultaneous normal stress and a shear stress will cause the crack to propagate at an angle 0 relative to the 
plane of the initial flaw or microcrack. The anisotropy of material properties can also affect the angle of 
propagation. This direction of growth will be the direction of cracking that would be observed in rupture 
experiments (or loading experiments to the point of reaching the proportional limit or some level of 
damage threshold). Actually observing the orientation of the critical flaw would be difficult. This report 
does not consider the topic of crack propagation in an anisotropic medium. Nonetheless, the methodology 
described in this report provides a mechanistic-based framework for which to describe orientation 
distribution for critical flaws. 


3.0 Examples 

3.1 Isotropic Material Critical Flaw Orientation Probability Density Distribution 

Fortran was used to implement a gaussian quadrature numerical integration algorithm to compute the 
generalized probability of failure equation (Eq. (33)) and the PDF equation for the orientation of critical 
flaws (Eq. (42)). This algorithm procedure is the same as that described in Nemeth et al. (2005), except 
that the integration is performed sequentially over smaller discrete sectors of the unit sphere. The numeri- 
cal output from the Fortran code was used with MATFAB (MathWorks) software to plot the PDD (which 
herein are the results from the PDF) over the unit sphere. Values for parameters related to critical strength 
anisotropy (Eqs. (29) and (30)) are sensitive to the number of sampling points (number of gauss points, 
NGP) used in the numerical integration. When a greater number of points are used, the numerical solution 
asymptotically approaches a single value. In the following plots of the unit sphere (Figs. 6 to 1 1), a 
smaller number of sampling points was used (NGP = 5) for computational efficiency. This somewhat 
affected the numerical values of parameter r T but not the trends seen in the figures. Results for Table 1 and 
Figures 12 to 15 were generated for NGP = 30, which represents fully converged solutions. 
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TABLE I.— VALUES OF PARAMETER r T ( SEE EQ. (30)) FOR VARIOUS STRENGTH RATIOS 
OF THE STRONGEST TO WEAKEST MATERIAL STRENGTH DIRECTIONS 
UNDER UNIAXIAL TENSILE LOADING 


Strength ratio 

Parameter r T 


Parameter, t, T , deg 


10 

90 1 


Shetty shear sensitivity, C 


1.1236 

100.0 

1.1236 

100.0 

1.25 

1.57 

1.47 

1.39 

1.29 

1.5 

1.95 

1.83 

1.82 

1.59 

2.0 

2.68 

2.52 

2.83 

2.21 

3.0 

4.15 

3.91 

5.44 

3.51 

4.0 

5.64 

5.32 

8.83 

4.87 

6.0 

8.68 

8.18 

17.90 

7.68 

8.0 

11.77 

11.10 

30.01 

10.59 

10.0 

14.91 

14.06 

45.18 

13.58 


Figure 6 shows an example of the PDD from Equation (42) of critical flaws on the unit sphere for an 
isotropic material with a 100-MPa uniaxial tensile load applied in various directions. Figures 6(a), (b), 
and (c) show results for a uniaxial load applied in the x, y, and z material axes, respectively. The plots 
show the direction normal to the flaw plane as depicted schematically in Figure 1. As would be expected, 
the “spot” is centered on the direction of applied load, with the highest probability of orientation 
perpendicular to the applied load (dark red (center of “bullseye” shapes) indicates the highest probability 
and blue (background color) indicates the lowest probability). Figures 6(a), (b), and (c) are shown to 
demonstrate that the numerical algorithm and plotting software correctly reproduce the same results 
regardless of the direction of the applied load. Figure 6(d) shows the PDD numerical distribution versus 
angles a and P for an applied uniaxial tensile load in the z direction. Figure 6 is generated for an 
incremental unit volume with a load level corresponding to P/y = 0.6321 (for c oV = 100 MPa-m 110 ), 
m i /= 10, and C = 1.5 for mild shear sensitivity (using Eq. (11)). All subsequent plots were generated for 
P f v= 0.6321 and o oV = 100 MPa-m 1/10 . __ 

The size (angular extent) of the “spot” (or distribution) shown in Figure 6 is a function of m v and C . 
This is demonstrated in Figure 7 for a uniaxial tensile load applied in the y direction. The plots show that 
the higher m v is, the more localized (more confined or highly aligned) the “spot” with the applied load is 
(flaws are more perpendicular to the applied load) as shown in Figures 7(a), (c), and (e). The same is also 
true with C , where for the shear-insensitive condition ( C = 100.0), the spot is more tightly distributed as 
shown in Figures 7(b), (d), and (f) and, as C gets smaller (shown as C =1.5 for mild shear sensitivity), 
the spot is less tightly distributed as shown in Figures 7(a), (c), and (e). 

Figure 8 shows the effect of shear-sensitivity parameter C for an applied shear stress x yz for in y = 10. 
Results for increasing shear sensitivity for C = 100.0 (shear-insensitive), C = 1.5 (mildly shear- 
sensitive), and 0.82 (highly shear-sensitive) are shown in Figures 8(a), (c), and (e), respectively. In 
Figures 8(a) and (b), only the normal tensile stress component contributes to failure, and the spot of 
highest probability is centered on a = 90° and P = 45°. In Figures 8(c) and (d), for the case of mild shear 
sensitivity, the spot of the highest probability bifurcates and the two spots are symmetrically offset at a 
small P angle from a = 90° and P = 45°. This indicates that the effect of shear sensitivity on the micro- 
crack alters the angle of most probable fracture for an applied shear stress. This is further shown in 
Figures 8(e) and (f) where a high shear sensitivity shows two distinctly different fracture planes, but their 
symmetry is oriented about a = 90° and P = 45°. An additional note about Figures 8(e) and (f) is that, 
when the normal stress on the flaw plane becomes compressive, the failure probability associated with 
that orientation for the tensile failure mode is set to zero in the numerical method; hence, there is some 
truncation that can be seen in the figure, where the color fringes drop out (left side of the left spot and 
right side of the right spot). 
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(c) -1.0 y (d) 

Figure 6. — Probability density distribution (PDD) plots of the orientation of critical flaws on the unit sphere 
for a 100-MPa uniaxial tensile load for an isotropic material. The center of the “spot” indicates the highest 
likelihood of orientation. This tests the consistency of calculated results for loading in various directions. 
Here the Weibull modulus, mv = 10.0, and the Shetty shear sensitivity, C =1.5. (a) Load in the x direc- 
tion. (b) Load in the y direction, (c) Load in the z direction, (d) PDD in terms of magnitude normalized to 
radians for angles a and p of the load in the z direction. 
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Figure 7. — Probability density distribution (PDD) plots of the orientation of critical flaws on the unit 
sphere for an isotropic material with a 100-MPa uniaxial tensile load applied in the y direction, 
showing the effect of Weibull modulus, mv, and Shetty shear sensitivity, C . (a) mv = 5.0 and 
C = 1.5. (b) m v = 5.0 and C =_100.0. (c) m v = 10.0 and C_= 1.5. (d )m v = 10.0 and 
C = 100.0. (e) mv - 20.0 and C =1.5. (f) m v = 20.0 and C = 100.0. 
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Figure 8. — Probability density distribution (PDD) plots of the orientation of critical flaws on the unit sphej^ for 
an isotropic material for a x yz = 100-MPa shear stress showing the effect of Shetty shear sensitivity, C . 

(a) Weibull modulus, m v = 10.Cf_and C = 100.0. (b) PDD in terms of magnitude for mv= 10.0 and 
C = 100.0. (c) m v = 10.0 and C = 1.5. (d) PDD in terms of magnitude for m v = 10.0 and C =1.5. 

(e) m v - 10.0 and C = 0.82. (f) PDD in terms of magnitude for m v - 10.0 and C =0.82. 
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Figure 9. — Probability density distribution (PDD) plots of the orientation of critical flaws on the unit sphere for 
an isotropic material, (a) Equibiaxial tensile stress applied in the y and z directions, a x = a y = 100 MPa. 

(b) PDD in terms of magnitude of equibiaxial tensile stress, (c) Uniaxial compressive stress applied in the 
x direction, a x = 500 MPa. (d) PDD in terms of magnitude of uniaxial compressive stress. For the tensile 
failure mode, (a) and (b), the Weibull modulus, m v = 10.0; a 0 v= 100 MPa (MPa.m 1/10 ), C = 1.5 (Eq. (11)), 
v = 0.22; and the probability of failure of material volume, Pw = 0.6321 . For the compressive failure mode, 

(c) and (d), using Equation (12), m eV = 10, the Weibull scale parameter for compression, a 0(? v= 500 MPa 
(MPa.m 1/10 ), and P N = 0.6321. 

Figures 9(a) and (b) show an example of the unit sphere PDD for critical flaw orientation for an 
equibiaxial tensile stress state applied in the y and z axial directions (cj x = 0; a y = cr z ). This plot was 
generated for a mildly shear sensitive material (C = 1.5). The critical flaws have equal likelihood to be at 
any angle p in the y-z plane. Figures 9(c) and (d) show an example of the unit sphere PDD for the critical 
flaw orientation for a uniaxial compressive stress state for the compressive stress failure criterion using 
Equation (12). The critical flaw planes are oriented at 45° to the direction of applied loading — the 
x direction in this case. 

3.2 Transversely Isotropic Material Critical Flaw Orientation Probability Density 
Distribution for Off-Axis Uniaxial Load 

In this section, some examples of the critical flaw orientation PDD are shown for the critical strength 
anisotropy model of Section 2.3.2. Results for the flaw orientation anisotropy model of Section 2.3.1 are 
not shown (because plotted results would appear similar to that of the critical strength anisotropy model). 
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Figure 10 is an example of the unit sphere PDD for critical flaw orientation for a 2: 1 transversely 
isotropic strength response. This model assumes critical strength anisotropy and uses Equation (30) to 
define the weaker directions of strength (the T distribution). In this case, a 50-MPa uniaxial tensile load 
(in the y material direction) is applied in the plane with the weakest material strength (the y-z plane) for 
Pjv= 0 .6321. Figures 10(a) and (b) show the unit sphere PDD for the critical flaw orientation for a shear- 
insensitive (C = 100.0) response, and Figures 10(c) and (d) show the results for a shear-sensitive response 
( C = 1.1236, where K lc = (K llc or K ]Ul; ) for v = 0.22), where the mode I strength of the flaw and the 
mode II (or mode III) of the flaw are equal. The strength response is for a material that has a more gradual 
critical strength transition from the strong to the weaker material strength direction (^ T = nil radians (90°) 
and y r = 1). As would be expected, a more oval or elliptical spot is observed in comparison to the 
isotropic material in Figure 6(b) because of the effect of the weaker y-z material plane. 



Figure 10. — Probability density distribution (PDD) plots of the orientation of critical flaws on the unit sphere for a 
shear-insensitive and shear-sensitive transversely isotropic material. This is for a 50-MPa uniaxial tensile load 
(in the y material direction) applied along the plane with the weakest material strength (the y-z plane). The strength 
ratio is 2:1, where critical mode I stress-intensity factor, K\ c , anisotropy functions representing the transverse 
(equatorial-belt) distribution, ^r = n/2 radians (90°) and yr = 1 . (a) Plot on unit sphere, (b) PDD in terms of magni- 
tude for Shetty shear sensitivity, C = 100.0 (shear-insensitive response), and parameter in K\ c anisotropy function, 
rr = 2.21 . (c) Plot on unit sphere, (d) PDD in terms of magnitude for C = 1.1236 (shear-sensitive response) and 
y t =2.83. In both cases, the Weibull modulus, mv= 10; the Weibull scale parameter, a 0 v = 1 00 MPa (MPa.m 1/10 ); 
and the probability of failure of material volume, Ptv - 0.6321 . 
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Figure 1 1 is an example of the unit sphere PDD for the critical flaw orientation for a 3: 1 transversely 
isotropic strength response for a uniaxial load offset 15° from the strong material direction (a = 15° and 
P = 0°). This model also uses Equation (30) to define the weaker directions of strength (the T distribu- 
tion); however, the strength response is for a material that has a more abrupt critical strength transition 
from the strong to the weaker material strength direction (<; r = 0.08727 radians (5°), y t = 0, and r r = 3.56). 
The load is a uniaxial tensile stress of 99.3 MPa for a P/v= 0 .6321. A relatively mild shear sensitivity 
( C = 1 .4) is assumed. The parameter y T = 0 corresponds to a uniform distribution over the equatorial 
belt where K lc is low relative to the rest of the unit sphere. This modeling scenario is arbitrary but may be 
considered to model the effect of a weak interface, such as between a fiber and matrix or within the 
matrix constituent, or possibly represent anisotropy in thin-film coatings where interfacial layers and 
anisotropically oriented microstructures are present. In addition to the PDD plot on the unit sphere 
(Fig. 1 1(a)), a numerical plot of the PDD versus angle of orientation of the flaw normal is shown 
(Fig. 1 1(b)). Contrasting the unit-sphere plots of Figures 10(a) and (c) with Figure 1 1(a), shows the 
effect of parameter £, T as a narrower strip (the smaller extent of angle a), where the critical flaw has a 
high likelihood of being centered on a = 90° and P = 0°. Figure 1 1 also shows the relative sensitivity 
of strength anisotropy with critical flaw angle. When the offset uniaxial load is only 15°, the probability 
that the critical flaw will be oriented about a = 90° and P = 0° becomes significant. This leads to the 
question of how anisotropy affects the orientation of the most critical flaw that is predicted with this 
model. 

Figures 12 to 15 were prepared to help investigate that question or to at least demonstrate that 
consequence with the unit-sphere model. These graphs plot the angle of highest probability (the 
orientation where the critical flaw PDD value is maximum) versus the angle of a uniaxial tensile load 
applied at some angle a. These figures were prepared using m v = 100 MPa-m 110 and Pjv= 0.6321. The 
strongest material direction is at a = 0°, and the weakest material direction is at a = 90°. It bears 



Figure 11. — Probability density distribution (PDD) of critical flaw normals for a 3:1 strength response for a trans- 
versely isotropic material. The load is a uniaxial tensile stress of 99.3 MPa applied at a 15° offset from the x axis 
(a = 15° and p = 0°). A relatively mild shear sensitivity ( C = 1.4) is assumed, (a) Plot on unit sphere, (b) PDD in 
terms of magnitude for Weibull modulus, mv= 10; Weibull scale parameter, a 0 v= 100 MPa.m 1/10 ; Shetty shear 
sensitivity, C =1.4; Poisson’s ratio, v = 0.22; parameters in the K\ c anisotropy function representing the transverse 
(equatorial-belt) distribution: ^r= 0.08727 radians (5°), y r = 0, and rj = 3.563; probability of failure of material 
volume, Pfv= 0.6321. 
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repeating that this is a prediction of the critically oriented flaw and not the direction of the subsequent 
crack propagation. The figures illustrate the effect of shear sensitivity C , the abruptness of the critical 
strength transition from strong to weak response using parameter £, r , and the strong direction to weak 
direction strength ratio modeled with parameter r T (see Table I). This line of inquiry was motivated by 
the observation that in CMCs the angle of microcracking for offset uniaxial tensile loads tends to be 
perpendicular to the direction of applied loading (e.g., Cady et al. (1995) and Lynch and Evans (1996)). 
However, this is not universally the case, as shown in Figure 16 for a 0°/90° CMC laminate with a 
30°-offset loading. In that case, the matrix cracks tended to be misaligned 12° to 15° with the perpen- 
dicular to the loading direction and approximately aligned symmetrically 45° from either fiber axis. 

Figure 16 shows that adjacent plies can influence the microcracking pattern. This suggests that it may 
be plausible to consider the influence of the laminate as a whole and to construct anisotropic unit-sphere 
matrix failure criteria appropriately. However, at the microscale level, CMC failure behavior involves 
more complicated mechanisms including fiber interface debonding and sliding resistances. Therefore, for 
a CMC, the anisotropic unit-sphere material failure criterion may prove useful at the macroscale of the 
composite material or alternatively at the microscale of the individual constituents involving the matrix 
with the fiber where multiscale modeling would be requisite. At the macroscopic scale, Figures 12 to 15 
would be appropriate for a unidirectional laminate where all the fibers are aligned in the same direction. 
When the individual plies are not aligned with one another, Equations (29) and (30) would have to be 
suitably modified. 



Figure 12. — Most probable orientation of the critical flaw relative 
to the offset angle of an applied uniaxial tensile load (at angle 
a on the unit sphere) for various anisotropic strength ratios. 
This example is for a mild critical strength transition (param- 
eters in the critical mode I stress-intensity factor, K tc , aniso- 
tropy function representing the equatorial-belt distribution, 

= n/2 radians (90°) and y T = 1 ) and high shear sensitivity 
C = 1.1236 (mode I and mode II strength are equal). The 
strongest material direction is at a = 0°, and the weakest 
material direction is at a = 90°. PDD, probability density 
distribution. 
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Figure 13. — Most probable orientation of the critical flaw relative 
to the offset angle of an applied uniaxial tensile load (at 
angle a on the unit sphere) for various anisotropic strength 
ratios. This example is for a mild critical strength transition 
(parameters in the critical mode I stress-intensity factor, K| C , 
anisotropy function representing the equatorial-belt distri- 
bution, ^r= rt/2 radians (90°) and yr = 1) and shear insensi- 
tivity C = 100.0. The strongest material direction is at a = 0°, 
and the weakest material direction is at a = 90°. PDD, 
probability density distribution. 


Figure 12 shows the most probable orientation of the critical flaw relative to the offset angle for an 
applied uniaxial load for a mild mode I critical strength transition (1- T = tc/2 radians (90°) and y T = 1) and 
high shear sensitivity C = 1.1236 (where mode I and mode II strength are equal). Figure 13 shows the 
results for the case of shear insensitivity C = 100.0 (with all other parameters the same). The mild critical 
strength transition is representative of materials such as graphite. Figures 12 and 13 indicate that the 
critical flaw angle deviates from the angle of offset load depending on the level of shear sensitivity and 
the degree of strength anisotropy. Figure 12 indicates that, when the uniaxial tensile load is applied in the 
strong material direction (a = 0°), the critical flaw is already oriented at a significant angle of a away 
from the direction of applied load. This contrasts with Figure 13, where this angle is smaller — particularly 
for lower values of strength ratio, and tends to zero at a = 0°. Flowever, this angle of deviation increases 
rapidly as the angle of offset load increases from a = 0° for the higher levels of strength ratio. Figure 13 
shows that the gradual critical strength transition from the strong material direction to the weak material 
direction strongly influences the orientation of the critical flaw for an offset load, and Figure 12 shows the 
additional influence of shear sensitivity combined with the mode I critical strength anisotropy. 
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Figure 14. — Most probable orientation of the critical flaw relative 
to the offset angle of an applied uniaxial tensile load (at angle 
a on the unit sphere) for various anisotropic strength ratios. 
This example is for a strong critical strength transition (param- 
eters in the critical mode I stress-intensity factor, K\ c , aniso- 
tropy function representing the equatorial-belt distribution, 

^r= Chi 745 radians (10°) and yr= 1) and high shear sensi- 
tivity C = 1.1236 (mode I and mode II strength are equal). 
The strongest material direction is at a = 0°, and the weakest 
material direction is at a = 90°. PDD, probability density 
distribution. 


Figure 14 shows the most probable orientation of the critical flaw relative to the offset angle for an 
applied uniaxial load for an abrupt mode I critical strength transition (£,r = 0.1745 radians (10°) and 
y T = 1) and high shear sensitivity, C = 1.1236 (where mode I and mode II strength are equal). Figure 15 
shows the results for the case of shear insensitivity C = 100.0 (with all other parameters the same). The 
abrupt mode I critical strength transition is intended to be representative of materials where an interface or 
inhomogeneity or abrupt shift in microstructure is present such as an interfacial boundary between the 
matrix and a fiber coating or the fiber itself, or between different deposition methods and/or materials on 
thin-film brittle coatings. 

Figures 14 and 15 show the abrupt transition of the orientation of the critical flaw to the weak 
material direction a = 90° as the angle of the offset load increases. The abruptness of this transition can 
also be seen in Figure 11, where the offset angle was chosen near to where this transition occurs. Both 
figures indicate expected isotropic behavior prior to this transition where the most probable angle of the 
critical flaw is coincident with the angle of the applied loading (as also seen in Fig. 6 for the different 
directions of applied loading). In Figure 14, the shear sensitivity shifts the angle of where this transition 
occurs to smaller values of a in comparison to the shear insensitive condition of Figure 15. As the ratio of 
the strength anisotropy increases, it also shifts this transition to progressively smaller values of a. 
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A comparison of Figures 12 and 13 with Figures 14 and 15 shows the effect of parameter 't T on the 
most probable angle of the critical flaw. Even for a milder transition of the mode I critical strength 
(£,r = tc/ 2 radians (90°)), the most probable orientation for the critical flaw trends toward a = 90° as the 
strength ratio increases. Overall, Figures 12 to 15 show the strong trend that the most probable orientation 
for the critical flaw is to shift toward the weaker material direction (a = 90°) as the strength anisotropy 
ratio increases. 



Figure 15. — Most probable orientation of the critical flaw relative 
to the offset angle of an applied uniaxial tensile load (at angle 
a on the unit sphere) for various anisotropic strength ratios. 
This example is for a strong critical strength transition (param- 
eters in the critical mode I stress-intensity factor, K\ c , aniso- 
tropy function representing the equatorial-belt distribution, 

^r= 0.1745 radians (10°) and yr= 1) and shear insensitivity 
C = 100.0. The strongest material direction is at a = 0°, and 
the weakest material direction is at a = 90°. PDD, probability 
density distribution. 


NASA/TM- 


-2013-217810 


29 



Figure 16. — Matrix cracking pattern from 30° off-axis loading of 
a 0°/90° ceramic matrix composite laminate from Lynch and 
Evans (1996). Matrix cracks appear to be 12° to 15° mis- 
aligned with the perpendicular to the loading direction, where 
a is the direction of loading, b is perpendicular to the loading 
direction, c is the fiber axis of the exposed ply, d is the fiber 
axis of the underlying ply, and e is 45° to either fiber axis. 


4.0 Summary and Conclusions 

A generalized equation was developed to describe the probability density function for the distribution 
of critical flaws from the Batdorf unit-sphere model for multiaxial loading. This formulation is for iso- 
tropic and anisotropic (transversely isotropic) brittle materials. The scenarios considered for anisotropic 
strength response were for (1) flaw orientation anisotropy, whereby a preexisting microcrack has a higher 
likelihood of being oriented in one direction than another direction, and (2) critical strength or fracture 
toughness anisotropy, where the level of critical strength or fracture toughness for mode I crack propaga- 
tion changes with regard to the orientation of the microstructure. 

A numerical algorithm was developed, and the results were plotted in MATLAB (MathWorks) for the 
unit-sphere probability density distribution (PDD) for various multiaxial stress states for an isotropic 
material. The results demonstrated that the PDD depends on the stress state, shear sensitivity of the flaws, 
and the Weibull modulus. Lower Weibull modulus and higher shear sensitivity meant higher dispersion of 
the distribution of critical flaws for a uniaxial stress state. For a pure shear stress state, when the shear 
sensitivity of the flaws increased, the orientation distribution of the critical flaws bifurcated into the two 
preferred angles of highest probability. This report also showed PDD results for an equibiaxial stress state 
and a uniaxial compressive stress state. 
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For transversely isotropic strength response, an example was shown of the PDD for a uniaxial tensile 
load applied at an angle offset from the strongest material direction. The angle of highest probability for 
the critical flaw also was investigated versus the offset loading angle for various anisotropic strength 
ratios and abruptness of the strength dropoff for the mode 1 critical strength anisotropy model. These 
results indicate that when there was an abrupt dropoff in mode 1 critical strength, there was also an abrupt 
transition of the angle of the critical flaw versus the offset angle of the applied uniaxial load such that the 
critical flaw angle became normal to the weakest material direction. These results also were sensitive to 
the level of shear sensitivity. Even for a more gradual dropoff in mode 1 critical strength (versus 
orientation) the most probable critical flaw angle tended to align with the weakest material direction as 
the loading offset angle increased and as the anisotropic strength ratio increased. Overall, there was a 
strong tendency for the most probable angle of the critical flaw to quickly align with the weakest material 
direction as the offset loading angle increased and as anisotropic strength ratio increased, even with a 
more gradual transition of mode 1 critical strength anisotropy versus orientation on the unit sphere. 

If one knows the orientation distribution of critical flaws (and the subsequent direction of crack 
propagation), the anisotropic stiffness degradation (the anisotropic elastic constants associated with the 
damaged material) can be determined. The author anticipates that this will be helpful for a follow-on 
phase of this effort, not described here, of enabling the unit-sphere failure criterion methodology to work 
with NASA’s micromechanics analysis code/generalized method of cells (MAC/GMC, Bednarcyk and 
Arnold (2002)). This incoiporation will allow the full exercise of the unit-sphere methodology, including 
incremental time/load steps and fatigue analysis (as described in Nemeth, Noel N.; Gyekenyesi, J.P.; and 
Jadaan, Osama M. (2005): Lifetime Reliability Prediction of Ceramic Structures Under Transient 
Thermomechanical Loads. NASA/TP — 2005-212505. http://ntrs.nasa.gov/), to predict the durability 
(strength and lifetime) of composite laminates and woven composite structures. 

The development of the anisotropic unit-sphere methodology was an attempt to provide an improved 
mechanistic basis to the problem of predicting composite strength under multiaxial loading in comparison 
to the phenomenologically based polynomial formulations such as Tsai-Wu, Tsai-Hill, and Hashin. The 
interesting consequence of this work was that the PDD on the unit sphere of the orientation of critical- 
strength- controlling flaws could also be predicted. That formulation and the demonstration of that 
capability was the subject of this report. 

This technology has general relevance to the anisotropic strength response of brittle materials. This 
may include thin- film coatings such as environmental barrier coatings where tensile stresses play a 
significant role in coating integrity. The need for additional enhancements is anticipated, such as the 
determination of the angle of noncoplanar crack growth and the development of the anisotropic stiffness 
matrix for a damaged material. 
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Appendix — Symbols 


A 

dA 

a 

a c 

C 

c 

m 

7L(wP) 

H 

K, 

K lc 

ksv 

kBV 

k w y 

m eV 

my 

n 

P(X = x) 

Pf 

Pjv 

A Pjy 

C Pjv)i 

Ps 

P sV 

(P sV )i 

XPxv 

Piv 

p 


area 

infinitesimal area on the surface of a unit radius sphere 

crack length or crack radius 

critical crack length or crack radius 

Shetty shear-sensitivity coefficient 

constant 

normalizing constant 

discrete density function of random variable X where P(X = x) 
normalized anisotropy function of K lc or as a function of angles a and P 

Heaviside step function 
mode 1 stress-intensity factor 
critical mode 1 stress-intensity factor 
Batdorf crack density coefficient 

Batdorf uniaxial stress-state normality constant for the volume-flaw 
failure mode (also known as the Batdorf normalized crack-density coefficient) 

uniaxial Weibull crack-density coefficient 

Weibull modulus for the volume-flaw failure mode for compressive stress states 
Weibull modulus for the volume-flaw failure mode 
number of incremental links; number of failure modes 

probability that a discrete real-valued random variable X equals a possible value x 
probability of failure (P/= 1 —P s ) 
probability of failure of material volume 

probability of failure of a crack with a strength between ai, and (gi, + Agi, ) in AK 

probability of failure of the z' th link 

reliability or probability of survival (P s = 1— Pf) 

probability of survival of material volume V 

probability of survival of the i link 

probability of existence of a crack with strength between Gi c and (G[, + Ag^ ) in an 
incremental volume 

probability that a crack of critical strength will be oriented in a particular direction 
such that it will grow and cause failure 

success probability 
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P(a, P) 

Pc( 2, a, P) 

Pz(a, P) 
£>z(iX, a, P) 

n 

r T 

s 

Y 

AV 
A Vi 
X 
x 

x, y, z 

Y 

a, P 

Y 
Yi 
It 

m 

C(a, P) 

lit' 

r|[{a, P) 

0 

A 


probability that a crack of critical stress <3 ieqc is oriented in the range between a and 
(a + da) and between p and (P + dp) 

probability that the critical flaw (the failure-initiating flaw) is oriented in the range 
between a and (a + da) and between p and (P + dp) under stress state X 

probability density distribution (PDD) function of fp (a, P) 

PDD function of the orientation of critical flaw normals from the applied stress state 
X given by angles a and P 

constant (ratio) or parameter in K lc anisotropy function for polar-cap distribution 
constant (ratio) or parameter in K lc anisotropy function for equatorial-belt distribution 
slope of a line 
volume 

incremental volume 

incremental volume associated with the i th li nk 
discrete real-valued random variable 

real valued number; possible value of discrete real-valued random variable X 
location in the body of the structure; Cartesian coordinates 
crack-shape geometry factor 

orientation angles or angular coordinates, where the direction normal to the plane of 
the microcrack is specified by the radial line defined by a and P in stress space 

constant or parameter in K ic anisotropy function representing the exponent of the sine 
or cosine function 

constant or parameter in K lc anisotropy function representing the longitudinal (polar- 
cap) distribution 

constant or parameter in K lc anisotropy function representing the transverse 
(equatorial-belt) distribution 

^ as a function of angle a describing the anisotropy of flaw orientation 

function describing the anisotropy of the flaw orientation where the normal direction 
to the flaw plane is given by angles a and P 

crack-density function 

crack-density function dependent on angles a and P 

crack-density function for applied uniaxial stress, where a is the number of flaws per 
unit volume with strength equal to or less than a 

crack-density function for equivalent mode 1 strength, G\ eqc , of a flaw: number of 
flaws per unit volume with strength equal to or less than a leqc 

angle of crack propagation relative to the initial plane of the microcrack 

constant or parameter in flaw-orientation anisotropy function representing one-half of 
the total angular extent of the anisotropy distribution 
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At constant or parameter in flaw-orientation anisotropy function representing the 

longitudinal (polar-cap) distribution 

A T constant or parameter in flaw-orientation anisotropy function or in K\ c anisotropy 

function representing the transverse (equatorial-belt) distribution 

A, positive number 

v Poisson’s ratio 

'c, L constant or parameter in K u anisotropy function representing the longitudinal (polar- 

cap) distribution 

Z, T constant or parameter in K lc anisotropy function representing the transverse 

(equatorial-belt) distribution 

I summation function; applied far-field multiaxial stress state 

g applied uniaxial stress 

G| r critical mode 1 strength 

a,,, x , CT, r y , o Ic z orthogonal critical strength components normalized by Gi c , max 


<3] eq equivalent, or effective, stress 

Gieg, max maximum value of G\ eq over the unit sphere from the applied multiaxial stress I 

G| m/ (x, y, z, a, P) equivalent, or effective, stress as a function of location x, y, z and orientation a, P 

Gi^ critical equivalent, or effective, stress 

c?i«,c(x, y, z, a, P) mode 1 far-field strength of a flaw located at coordinates x, y, and z and oriented at 
angles a and P 

g, local uniaxial stress 

g„ applied far-field stress component normal to a crack face 

G„ Weibull scale parameter 

a oe v Weibull scale parameter for compression for the volume-flaw failure mode 

normalized to unit volume 

G 0 j/ Weibull scale parameter for the volume-flaw failure mode normalized to unit volume 

G„;/ threshold stress parameter 

g x G y , g z orthogonal stress components expressed relative to a global coordinate system 

t shear stress acting on the oblique plane whose normal is determined by angles a and 

P, which represents the applied far-field shear stress on a crack face 

T xy , Xy Z , Tzx shear stress components expressed relative to a global coordinate system 

<j) constant or parameter in flaw-orientation anisotropy function representing the 

exponent of the sine or cosine function 

<j>i constant or parameter in flaw-orientation anisotropy function representing the 

longitudinal (polar-cap) distribution 
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(j) r constant or parameter in flaw-orientation anisotropy function representing the 

transverse (equatorial-belt) distribution 

Q(X, G\ rqc ) area of a solid angle projected onto a unit radius sphere in three-dimensional stress 

space for which a lnj > <3i eqc from an applied multiaxial stress state I 


Superscript 

normalized 


Subscripts 
I, II, III 

i 

L 

max 

n 

T 

V 


mode I, II, or III 
i th value or I th term 
longitudinal 
maximum 

integer; number of links; normal 
transverse 

volume or a volume-based property (e.g., indicates volume-flaw analysis) 


Definitions 

L 

mode I 
mode II 
mode III 
T 


longitudinal 
crack-opening mode 
crack-sliding mode (in-plane shear) 
crack-tearing mode (out-of-plane shear) 
transverse 
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